@asalmgren The new review issues should be addressed now.
The deviations from the existing zlo path are fixed and checks are added to ensure other setups (SHOC/PBL) use MOST only on the zlo surface. The new ABL_MOST_Cloudchamber test setup is intentionally unchanged to mirror a real setup (using regional averaging and explicit zref on each side). However, unit tests were added to check for both planar and regional averaging with default zref so the missing coverage should be picked up in CI now.

Summary

Things I'd want addressed before merge

1. Plane averaging double-counts on lateral faces under tiling.

compute_plane_averages, first loop (ERF_MOSTAverage.cpp:1213-1228): the face test is on vbx (the grid), then every tile of that grid gets pbx.setSmall(dir,sm_index); pbx.setBig(dir,sm_index). TileNoZ() returns (mfiter_tile_size[0], mfiter_tile_size[1], 1024000), and AMReX's default mfiter_tile_size is (1024000,8,8) — so on CPU there is tiling in y. For a ylo/yhi face, every y-tile of a grid collapses onto the same y-plane and gets summed, while ncell_plane counts the plane once. A grid 32 cells deep in y over-counts by 4×. x-faces are safe only because the default x tile size is huge; set fabarray.mfiter_tile_size and they break too. The Tv block and the Umag block right below both test pbx instead and are correct — that internal inconsistency is the tell. Switching the first loop to the pbx test matches them and fixes it.

The new gtest wouldn't catch this: the domain is (3,5,7)–(5,8,11), one box, 4 cells in y, so no tile split. Bumping that domain past 8 cells in y (or adding a multi-box BoxArray) would cover it, and would also give the first multi-grid coverage of these paths.

2. t31_arr / t32_arr null-deref guard is inverted. ERF_SurfaceLayer.cpp:~690:

cpp
auto t31_arr = (dir == 0 || Tau_lev[TauType::tau31]) ? Tau_lev[TauType::tau31]->array(mfi) : Array4<Real>{};

When dir == 0 the short circuit skips the pointer check entirely, so a null tau31 dereferences. In practice update_diffusive_arrays allocates tau21/31/32 whenever l_Surf_X || l_Surf_Y, so this only fires with l_use_diff == false — but then it's a segfault where an assert would be much kinder. (Tau_lev[...] != nullptr) alone is the right condition; the dir == 0 clause buys nothing.

3. Ghost values on lateral faces are never communicated.

update_fluxes now skips FillBoundary on u_star/t_star/q_star/olen unless coordDir() == 2. But compute_u_flux for dir == 0 reads u_star_arr(ic,j-1,kc) and umm_arr(ic,j-1,kc) at the low y-node of every grid. With more than one grid in y along an x-wall, that's a ghost cell filled only by whatever compute_fluxes happened to compute over its own growntilebox, using average-MF ghosts that (see 4) are themselves not filled tangentially. The comment says the faces are "filled individually on corresponding ranks," which is true of the valid region but not the halo. Worth either restoring FillBoundary for all faces or documenting why the grown-tilebox computation is sufficient.

4. Periodicity(IntVect::TheDimensionVector(dir)) on the average FillBoundary calls (lines 1677, 1800, 1944-1947).

Periodicity takes period lengths, so this asks for period 1 in the collapsed direction and non-periodic in the other two. It does propagate the single valid plane into the z-ghosts, which I assume is the intent — but it also drops the tangential periodic images that geom.periodicity() used to supply. For a zlo run with periodic x/y and average_policy = 1/2, the ghost values feeding um_arr(i-1,...) at the periodic boundary now come from nowhere. No gold file changed, so either the affected configs aren't in CI or the ghosts don't reach the answer — but the reasoning should be in a comment rather than left implicit.

6. Cross-BoxArray MFIter indexing.

compute_fluxes and fill_qsurf_with_qsat now iterate *m_lmask_lev[lev][0] and call ->array(mfi) on MultiFabs built on a different BoxArray (collapsed in dir, not in z). This works only because both derive from the same 3D grids with the same DM and because of the new full-z-column assert — AMReX will not catch a violation, even in debug. The rationale in the comment is good; I'd add an AMREX_ASSERT on boxArray().size() and DM equality so a future decomposition change fails loudly. Relatedly, gtbx.grow(2,3) in fill_qsurf_with_qsat hardcodes a ghost width that must not exceed ng[0] of the 2D MF — worth deriving from nGrowVect().

@AMLattanzi @asalmgren all of the review comments should be addressed now

Review

Blocking

1. BULK_COEFF is reachable on a lateral face and is wrong there. In compute_SurfaceLayer_bcs (ERF_SurfaceLayer.cpp:1232-1264), the block uses klo = sm_index as a z index (cons_arr(i,j,klo), t13_arr(i,j,klo), hfx3_arr(i,j,klo)) and u_star_arr(i,j,0). On an x/y face sm_index is an x/y index, t13/hfx3 aren't written at all, and the 2D parameter fabs are collapsed in x/y rather than z — so erf.xlo.surface_layer.flux_type = bulk_coeff reads garbage and, for a high face, out of bounds. flux_type is parsed per face, so nothing currently prevents it. The PR description says only the surface-temperature + Moeng pathway is face-aware; that should be an AMREX_ALWAYS_ASSERT in the SurfaceLayer constructor (dir == 2 && isLow() for BULK_COEFF, and probably CUSTOM/RICO too), not just prose.

2. Asymmetry in the sampling of velocity at the high face. Consider set_k_indices_N, line 750:

const int lk = is_lo_face ? dom_lo + wall_offset : dom_hi - wall_offset;

With the default zref = dx/2, wall_offset = 0, so on xhi m_k_indx = dom_hi and on xlo m_k_indx = dom_lo. Since the u velocity lives at the face, the true sample point on the high face is dom_hi + 1. Please examine the code for other asymmetries like the one identified here.